﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.NonlinearEquationSolvers
{
    /// <summary>
    /// This class provides methods to find the roots of a polynomial by using the Bairstow
    /// root finding algorithm. 
    /// </summary>
    [Serializable]
    public class BairstowComplexRootFinder
    {
        /// <summary>
        /// First working array for the calculation.
        /// </summary>
        private readonly double[] a;

        /// <summary>
        /// Second working array for the calculation.
        /// </summary>
        private readonly double[] b;

        /// <summary>
        /// Third working array for the calculation.
        /// </summary>
        private readonly double[] c;

        /// <summary>
        /// Describes if the calculated roots have a precision error.
        /// </summary>
        private bool precisionError;

        /// <summary>
        /// The polynomial for calculating the roots.
        /// </summary>
        private readonly SimplePolynomial polynomial;

        /// <summary>
        /// Initializes a new instance of the <see cref="BairstowComplexRootFinder"/> class.
        /// </summary>
        /// <param name="polynomial">The polynomial for calculating the roots. The minimum 
        /// degree of the polynomial has to be three.</param>
        public BairstowComplexRootFinder(SimplePolynomial polynomial)
        {
            if (polynomial == (SimplePolynomial) null)
            {
                throw new ArgumentNullException("polynomial");
            }

            if (polynomial.Degree <= 2)
            {
                throw new ArgumentException("The minimum degree of the polynomial has to be three.");
            }

            this.a = new double[polynomial.Degree + 1];
            this.b = new double[polynomial.Degree + 1];
            this.c = new double[polynomial.Degree + 1];
            this.polynomial = polynomial.Copy();
        }

        /// <summary>
        /// Describes if the calculated roots have a precision error. If true the
        /// calculate roots may not be correct.
        /// </summary>
        /// <value><c>True</c> if [precision error]; otherwise, <c>false</c>.</value>
        public bool PrecisionError
        {
            get { return precisionError; }
        }

        /// <summary>
        /// Finds the roots of the polynomial.
        /// </summary>
        /// <returns>The complex roots of the specified polynomial.</returns>
        public ComplexNumber[] FindRoots()
        {
            int n = this.polynomial.Degree;
            List<ComplexNumber> result = new List<ComplexNumber>();

            this.b[0] = 1;
            this.c[0] = 1;
            this.precisionError = false;
            this.polynomial.Normalize();

            for (int i = this.polynomial.Degree; i >= 0; i--)
            {
                this.a[i] = this.polynomial.ValueAtRank(this.polynomial.Degree - i);
            }

            while (n > 2)
            {
                this.InternalRootFinder(n);

                ComplexNumber[] res = ExtendedMath.SolveQuadraticEquationComplex(this.a[n - 1], this.a[n]);

                result.Add(res[0]);
                result.Add(res[1]);

                n -= 2;
            }

            if ((n % 2) == 1)
            {
                result.Add(new ComplexNumber(-this.a[1]));
            }
            else
            {
                ComplexNumber[] res = ExtendedMath.SolveQuadraticEquationComplex(this.a[n - 1], this.a[n]);

                result.Add(res[0]);
                result.Add(res[1]);
            }

            return result.ToArray();
        }

        /// <summary>
        /// The internal algorithm to find the roots.
        /// </summary>
        /// <param name="n">The nth root to find.</param>
        private void InternalRootFinder(int n)
        {
            double r = 0;
            double s = 0;
            double dr = 1.0;
            double ds = 0;
            double epsilon = 1e-14;
            int iterations = 1;
            Random rand = new Random();

            while ((ExtendedMath.Fabs(dr) + ExtendedMath.Fabs(ds)) > epsilon)
            {
                if ((iterations % 200) == 0)
                {
                    r = (double) rand.Next(32767) / 16000;
                }

                if ((iterations % 500) == 0)
                {
                    epsilon *= 10.0;
                    precisionError = true;
                }

                b[1] = a[1] - r;
                c[1] = b[1] - r;

                for (int i = 2; i <= n; i++)
                {
                    b[i] = a[i] - r * b[i - 1] - s * b[i - 2];
                    c[i] = b[i] - r * c[i - 1] - s * c[i - 2];
                }

                double dn = c[n - 1] * c[n - 3] - c[n - 2] * c[n - 2];
                double drn = b[n] * c[n - 3] - b[n - 1] * c[n - 2];
                double dsn = b[n - 1] * c[n - 1] - b[n] * c[n - 2];

                if (ExtendedMath.Fabs(dn) < 1e-16)
                {
                    dn = 1;
                    drn = 1;
                    dsn = 1;
                }

                dr = drn / dn;
                ds = dsn / dn;
                r += dr;
                s += ds;
                iterations++;
            }

            for (int i = 0; i < n - 1; i++)
            {
                a[i] = b[i];
            }

            a[n] = s;
            a[n - 1] = r;
        }
    }
}